Naval  Research  Laboratory 

Washington,  DC  20375-5320 


NRL/MR/6794-97-7893 


Electrostatic  Particle-in-Cell 
Simulation  Technique  for 
Quasineutral  Plasma 

Glenn  Joyce 
Martin  Lampe 
Steven  P.  Slinker 
Wallace  M.  Manheimer 

Beam  Physics  Branch 
Plasma  Physics  Division 


March  31,  1997 


Approved  for  public  release;  distribution  unlimited. 


19970410  043 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  Information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
collection  of  Information  includino  supqestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson 


1.  AGENCY  USE  ONLY  {Leave  Blank) 

2.  REPORT  DATE 

March  31,  1997 

3.  REPORT  TYPE  AND  DATES  COVERED 


Interim  Report 


4.  TITLE  AND  SUBTITLE 


Electrostatic  Particle-in-Cell  Simulation  Technique  for  Quasineutral  Plasma 


6,  AUTHOR(S) 

Glenn  Joyce,  Martin  Lampe,  Steven  P.  Slinker,  and  Wallace  M.  Manheimer 


7.  PERFORMING  ORGANIZATION  NAME{S)  AND  ADDRESS(ES) 

Naval  Research  Laboratory 
Washington,  DC  20375-5320 


9.  SPONSORING/MONITORING  AGENCY  NAME{S)  AND  ADDRESS(ES) 
Office  of  Naval  Research 

Arlington,  VA  22217 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

NRL/MRy6794-97-7893 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 


12a.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited. 


1 12b.  DISTRIBUTION  CODE 


13.  ABSTRACT  {Maximum  200  words) 

We  have  developed  a  particle-in-cell  simulation  method  in  which  the  electrostatic  field  is  determined  from  the  requirement 
of  quasineutrality,  rather  than  Poisson’s  equation.  Time  steps  may  be  orders  of  magnitude  longer  than  the  plasma  period,  and 
mesh  cells  orders  of  magnitude  longer  than  the  Debye  length,  since  electron  plasma  oscillations  do  not  appear  in  the  model  and 
the  Debye  length  is  essentially  set  to  zero.  The  quasineutral  approach  also  avoids  the  problem  of  statistical  fluctuations  in  the 
charge  density,  which  frustrate  the  use  of  Poisson’s  equation  in  a  quasineutral  plasma.  The  simulation  technique  correctly 
represents  kinetic  features  such  as  Landau  damping.  The  method  is  demonstrated  by  application  to  several  simple  test  problems, 
including  free  expansion  of  a  plasma,  and  linear  and  nonlinear  ion  sound.  In  the  case  of  a  plasma  with  strongly  magnetized 
electrons,  we  apply  the  technique  to  determine  the  parallel  electric  field  and  parallel  transport  within  the  plasma.  Quasineutral 
techniques  for  representing  cross-field  transport,  and  edge  effects  in  bounded  plasmas,  will  be  discussed  in  subsequent 
publications. 


14.  SUBJECT  TERMS 

Quasineutral  plasma 

Plasma  simulation 

Particle  simulation 

Particle-in-cell  simulation 

17.  SECURITY  CLASSIFICATION 

OF  REPORT 

18.  SECURITY  CLASSIFICATION 

OF  THIS  PAGE 

19.  SECURITY  CLASSIFICATION 

OF  ABSTRACT 

UNCLASSIFIED 

UNCLASSIFIED 

UNCLASSIFIED 

NSN  7540-01-280-5500 


15.  NUMBER  OF  PAGES 

38 _ _ 

16.  PRICE  CODE 

20.  LIMITATION  OF  ABSTRACT 


Standard  Form  298  (Rev.  2-89) 
Prescribed  by  ANSI  Std  239-1 8 
298-102 


CONTENTS 


1.  INTRODUCTION .  1 

2.  CALCULATION  OF  THE  ELECTRIC  FIELD  .  3 

A.  Unmagnetized  One-Dimensional  Plasma  .  3 

B.  Strongly  Magnetized  Electrons .  6 

3.  FORMAL  ANALYSIS  OF  MODE  STRUCTURE  AND  STABILITY  .  7 

A.  Linearized  Normal  Modes .  ^ 

B.  Effect  of  Spatial  Smoothing .  O 

4.  NUMERICAL  IMPLEMENTATION .  13 

5.  AN  EXAMPLE:  FREE  EXPANSION  OF  A  PLASMA  .  15 

A.  Analytic  Treatment:  Exact  Quasineutrality  .  15 

B.  Analytic  Treatment  Using  Our  Model  .  16 

C.  Numerical  Simulation  .  1^ 

6.  SECOND  EXAMPLE:  ION  SOUND .  18 

A.  Standing-Wave  Ion  Sound  in  the  Linear  Regime .  19 

B.  Nonlinear  Traveling  Wave  .  20 

7.  CONCLUSIONS  .  21 

Appendix — VLASOV  THEORY  OF  ION  SOUND .  23 

REFERENCES  .  27 


iii 


ELECTROSTATIC  PARTICLE-IN-CELL  SIMULATION  TECHNIQUE 

FOR  QUASINEUTRAL  PLASMA 

1.  Introduction 

In  particle-in-cell  (PIC)  simulations  of  plasmas,  the  standard  technique*'^  for 
calculating  the  electrostatic  field  is  to  solve  Poisson’s  equation,  with  the  charge  density 
source  term  determined  by  the  laydown  of  the  densities  of  electrons  and  ions,  ne  and  n;. 
This  procedure  works  well  when  the  phenomena  of  interest  proceed  on  time  scales 
comparable  to  the  electron  plasma  frequency  coj,  and  spatial  scales  comparable  to  the 
Debye  length  JId.  and  when  there  is  substantial  charge  separation  between  electrons  and 
ions.  However,  there  are  many  plasma  problems  where  the  time  and  space  scales  are  very 
much  longer,  and  where  the  plasma  maintains  quasineutrality  throughout,  i.e.  Ine-nil «  Uj. 
In  these  situations,  it  can  be  extraordinarily  inefficient,  or  even  unfeasible,  to  use 
Poisson’s  equation. 

One  problem  is  that  the  simulation  supports  electron  plasma  oscillations,  and 
therefore  the  time  step  must  be  less  than  2(0j,'‘  for  stability.'^  In  many  situations,  electron 
oscillations  play  no  role  in  the  phenomena  of  interest,  and  the  shortest  time  scale  that  is 
actually  of  interest  may  be  the  period  of  a  low-frequency  wave,  an  ion  time  scale,  a 
collisional  time  scale,  or  a  time  scale  for  electron  transport  over  some  macroscopic 
length.  These  time  scales  may  be  several  orders  of  magnitude  longer.  Implicit 
algorithms^**  have  often  been  used  to  avoid  resolution  of  the  electron  plasma  oscillation 
time  scale.  However,  PIC  simulation  techniques  also  face  a  more  fundamental  difficulty 
arising  from  the  circumstances  of  quasineutrality.  The  charge  separation  between 
electrons  and  ions,  n*-  ni,  is  often  less  than  10'^  of  the  density  of  either  electrons  or  ions. 
If  the  electrons  and  ions  are  represented  by  simulation  macroparticles,  any  attempt  to 
calculate  the  potential  directly  from  Poisson’s  equation  would  be  futile,  and  overwhelmed 
by  statistical  noise.  For  example,  in  a  million-particle  2D  simulation  with  a  100  X  1(X) 
grid,  there  are  typically  100  macroparticle  electrons  or  ions  in  each  cell.  The  physically 
correct  value  of  the  difference  between  the  number  of  electrons  and  ions  in  the  cell  would 
be  on  the  order  of  10'^  macroparticles,  clearly  unresolvable,  whereas  the  statistical 

fluctuations  within  the  cell  would  be  on  the  order  of  VlOO  macroparticles,  which  is  four 
orders  of  magnitude  larger  than  the  actual  value.  Numerical  schemes  involving  Poisson’s 
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equation  are  obviously  very  difficult  (and  actually  inappropriate)  in  the  quasineutral  limit. 
Indeed,  Chen‘S  noted  long  ago  that,  “In  a  plasma,  it  is  usually  possible  to  assume  iie  =  ni 
and  VE^iO  at  the  same  time.  This  is  a  fundamental  trait  of  plasmas,  one  which  is 
difficult  for  the  novice  to  understand.  Do  not  use  Poisson’s  equation  to  obtain  E  unless  it 
is  unavoidable!” 

Over  the  years,  this  advice  has  been  applied  in  many  analytic  and  numerical 
models  which  represent  the  plasma  as  a  fluid,  or  represent  the  electrons  as  either  a 
dielectric  medium  or  a  fluid  within  some  hybrid  scheme.'^  *^  Methods  have  been 
developed  which  circumvent  the  use  of  Poisson’s  equation  by  neglecting  electron  inertia 
and  determining  E  from  the  resulting  simplified  electron  momentum  conservation 
equation.  In  these  models,  ni  is  determmed  from  d3mamical  equations,  but  iic  is  not 
separately  calculated;  it  is  simply  set  equal  to  ni  to  maintain  quasineutrality.  This 
procedure  eliminates  temporal  scales  on  the  order  of  the  electron  plasma  frequency,  as 
well  as  spatial  structures  on  the  Debye  length  scale.  However,  in  fully  kinetic  models, 
and  in  particular  in  PIC  simulations,  the  electron  density  is  calculated  directly  by  the 
simulation,  and  the  usual  approach  has  been  to  calculate  ^  from  Poisson’s  equation. 

We  have  developed  a  new  approach  to  the  simulation  of  quasineutral  plasmas  with 
particle  electrons  and  particle  ions.  Our  approach  is  motivated  by  the  quasineutral  fluid 
techniques  described  in  the  previous  paragraph,  and  our  objective  is  to  use  grid  spacing 
wide  compared  to  the  Debye  length,  and  time  steps  long  compared  to  the  electron  plasma 
frequency.  We  believe  that  the  technique  can  be  used  to  treat  a  wide  variety  of  plasma 
problems.  However,  the  primary  objective  of  our  present  work  is  multi-dimensional 
overall  modeling  of  an  electron  cyclotron  resonance  (ECR)  reactor  used  for  plasma 
processing.  The  plasma  in  this  case  is  bounded,  partially  ionized,  collisional,  and 
magnetized,  with  a  typical  plasma  density  lO'^  cm'^,  electron  temperature  several  eV,  and 
neutral  density  several  times  lO'^.  The  scale  size  of  the  reactor  is  tens  of  cm  (>10^A.d). 
and  the  time  scales  of  interest  range  from  10  ns  (electron  transit  times  over  cm  size 
features,  and  electron  collision  times)  to  hundreds  of  ps  (chemical  equilibration),  while 
X,D  =  10'^  cm  '  =  10'"  sec.  In  recent  years,  there  has  been  considerable  interest  in  the 

•  16-25 

use  of  particle-in-cell/Monte  Carlo  (PIC/MC)  codes  to  model  this  type  of  plasma. 
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Like  pure  PIC  codes,  these  have  typically  used  Poisson’s  equation  to  determine  the 
electric  field.  In  the  present  paper,  we  discuss  only  the  method  for  determining  the 
internal  electric  field  within  an  unmagnetized  bulk  plasma,  or  parallel  to  the  magnetic 
field  in  a  magnetized  bulk  plasma.  In  subsequent  publications,  we  shall  discuss  self- 
consistent  techniques  for  dealing  with  sheaths,  collisions  and  chemistry,  and  with  cross¬ 
field  transport  in  a  magnetized  plasma,  within  a  multidimensional  quasineutral 
framework. 

2.  Calculation  of  the  Electric  Field 

A.  Unmagnetized  One-Dimensional  Plasma 

For  simplicity,  we  shall  consider  a  one-dimensional  system  specified  by  Cartesian 
coordinate  z,  although  the  formalism  can  be  extended  to  multi-dimensional  systems.  We 
assume  that  the  simulation  is  globally  quasineutral,  i.  e.  the  total  number  of  electrons  is 
equal  to  the  total  number  of  ions.  We  also  assume  that  the  electron  Debye  length  is  small 
compared  to  any  scale  length  resolved  in  the  model,  and  the  electron  plasma  fiiequency  is 
fast  compared  to  any  time  scale  resolved  in  the  model.  Thus,  if  there  were  any  departure 
from  local  quasineutrality,  the  resulting  strong  electric  field  would  drive  electron  currents 
to  restore  quasineutrality  within  a  time  scale  of  several  electron  plasma  periods,  i.e., 
essentially  instantaneously  on  the  time  scale  of  the  model.  Thus,  the  electric  field  always 
takes  the  value  necessary  to  keep  the  electron  density  n*  equal  to  the  ion  density  ni.  To 
specify  this  electric  field,  we  can  begin  with  the  electron  momentum  conservation 
equation, 


1  9Pe 

-eE= - ^ 


(1) 


Here,  Pe(z,t),  Ue(z,t)  and  Ve(z,t)  are  the  electron  kinetic  pressure  (including  flow  terms), 
the  electron  fluid  velocity,  and  the  mean  electron  momentum  transfer  collision  frequency, 
which  can  be  specified  as  integrals  over  the  electron  distribution  function  fe(z,v,t). 
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Pe(z,  t)  S  J dv rngV^feCz,  V,  t). 


(2a) 


Ue(z,t)  =  ldvvfe(z,V,t), 

(2b) 

Ve(z,  t)  s  Jdv  V(Z,  V,  t) fc(z,  V,  t)  . 

(2c) 

These  integrals,  which  appear  in  the  first  two  terms  on  the  right  hand  side  of  Eq.  (1),  can 
be  evaluated  at  each  point  of  the  grid,  by  laying  down  the  mean  quantities  for  the 
electrons  assigned  to  that  grid  point.  [The  collision  frequency  v(z,v,t)  represents  a  sum 
over  the  various  collisional  processes  which  are  represented  in  the  simulation  as  Monte 
Carlo  events,  dynamical  friction,  etc.] 

The  first  two  terms  on  the  right-hand  side  of  (1)  represent  the  ambipolar  electric 
field,  which  balances  the  electron  pressure  gradients,  flow  gradients,  and  frictional  forces. 
For  low-frequency  plasma  processes  which  maintain  quasineutrality  (the  only  type  of 
processes  we  wish  to  follow),  the  inertial  term  [last  term  of  Eq.  (1)]  is  smaller  by  order 
me/mi.  hi  a  quasineutral  formulation,  the  inertial  term  would  be  neglected,  and  the 
ambipolar  field  would  be  used  in  the  ion  dynamical  equations  to  calculate  ni  at  the  next 
time  step.  Then  iic  would  simply  be  set  equal  to  Ui.  However,  this  is  not  quite  sufficient 
in  a  particle  simulation,  where  the  electrons  and  ions  evolve  separately  during  each  time 
step.  With  E  set  equal  to  the  ambipolar  field,  nc(z,t)  remains  constant  in  time  (except  for 
statistical  fluctuations),  while  the  ion  density  ni(z,t)  gradually  evolves  because  of  the 
relatively  slow  ion  motion.  Thus  the  ambipolar  field  alone  will  not  maintain  the 
quasineutrality  relation 


ne(z,t)  =  ni(z,t).  (3) 

Even  worse,  particle  simulations  are  always  subject  to  statistical  fluctuations  in  the 

density  and  flux  of  any  species  at  any  given  grid  point,  typically  of  the  order  of  Vn  , 
where  N  is  the  number  of  particles  in  a  grid  site.  Thus,  a  particle  simulation  code  must 


4 


contain  some  mechanism  for  stably  maintaining  the  quasineutrality  relation  (3)  in  the  face 
of  these  fluctuations,  which  are  of  much  larger  order  than  me/mi. 

We  have  found  that  for  low-frequency  phenomena  the  kinetic  information  of 
interest  is  essentially  contained  in  the  ambipolar  field,  and  that  the  last  term  of  (1)  serves 
only  the  functional  purpose  of  keeping  ne  equal  to  ni.  This  opens  up  the  possibility  of 
keeping  the  ambipolar  field  but  using  an  approximate  technique  to  maintain  Eq.  (3), 
rather  than  calculating  the  actual  inertial  term.  In  earlier  work,^^  we  have  experimented 
with  the  technique  of  pushing  the  both  the  electrons  and  ions  in  the  ambipolar  field  alone, 
and  then  applying  an  approximate  correction  field  to  the  electrons  which  restores  the 
quasineutrality  condition.  This  worked  well,  and  we  were  able  to  prove  that  kinetic 
properties  such  as  the  Landau  damping  of  ion  sound  waves  were  preserved.  However, 
this  approach  can  become  somewhat  complicated,  especially  when  care  is  taken  to 
conserve  energy  exactly.  In  the  present  paper,  we  describe  an  approach  which  is  even 
simpler,  works  extremely  well  over  very  long  times,  and  has  excellent  stability  and  energy 
conservation  properties.  This  approach  is  simply  to  replace  Eq.  (1)  with  a  modified  form 
of  the  ambipolar  field, 

1  d 

-eE  =  — ^niTj  +  Vem^Ue,  (4) 

nj  dz'  ' 

where  the  electron  kinetic  pressure  Pe(z,t)  =  ne(z,t)Te(z,t)  is  replaced  by  ni(z,t)Te(z,t), 
using  the  ion  density  instead  of  n*.  Here,  Te  is  a  kinetic  temperature  (including  flow 
terms),  defined  as 

ne(z,  t)Te(z,  t)  S  J  dvm^v^feCz,  V,  t) .  (5) 

This  simple  artifice  causes  the  electron  density  to  remain  closely  coupled  to  the  ion 
density.  This  can  be  seen  by  writing  the  electron  momentum  conservation  equation  in  the 
form 
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(6) 


1  9 

+  — +  eE  +  =  0. 


3u, 


3t  rig  3z 


Using  Eq.  (4)  for  E,  this  gives 


3U  rr.  ^  /I  I 

m--:^  =  T  — £n 

®  at  ®  az 


(7) 


Equation  (7)  shows  that  the  electrons  are  always  accelerated  up  the  gradient  nj/ne  i.e. 
toward  the  point  of  maximum  positive  charge  density.  The  result  is  that  the  electron 
density  oscillates  about  the  ion  density.  Although  these  oscillations  are  unphysical,  they 
are  rapid  and  the  oscillation  amplitude  should  remain  small  if  the  system  is  started  in  a 
quasineutral  state  and  the  time  step  is  small  enough.  The  stability  properties  of  these 
oscillations  will  be  discussed  in  Sec.  3,  and  examples  will  be  given  in  Secs.  5  and  6. 

B.  Strongly  Magnetized  Electrons 

In  the  application  which  we  are  studying,  ECR  plasma  sources,  the  electrons  are 
strongly  magnetized,  with  gyrofrequencies  comparable  to  coj;,  and  gyroradii  comparable  to 
the  Debye  length.  Since  we  do  not  wish  to  resolve  these  short  time  and  space  scales,  it  is 
convenient  to  use  a  guiding  center  representation  of  the  electrons.  An  electron  is 
characterized  by  its  coordinate  the  curvilinear  coordinate  along  the  field  line,  its 
parallel  velocity  vn  s  d^dt,  and  the  magnitude  vj.  of  its  perpendicular  velocity.  However, 
in  practice  is  more  convenient  to  use  the  electron’s  magnetic  moment  |i  =  mvx^/2B  as  the 
independent  variable,  rather  than  vx,  since  p  is  an  adiabatic  invariant.  Here,  B  is  the 
magnitude  of  the  magnetic  field.  The  equation  of  motion  for  an  electron,  between 
collisions,  includes  a  mirror  force  term  and  is 


dV|i 

dt 


=  -eE„-p-^. 


(8) 


6 


Within  our  quasineutral  formulation,  the  electric  field  component  parallel  to  B  is  then 
specified  by  the  electron  momentum  equation  in  the  form 


B  a  Pell  -dB 


+  Verne 


1  a  /  \ 

Uell+— ^KmeUellj. 


(9) 


where  Pen  and  Pex  are  the  electron  kinetic  pressure  parallel  and  perpendicular  to  B,  and  \l 
=  Tei/B  is  the  mean  magnetic  moment  for  electrons  at  a  given  location,  [Note  that  the 
mirror  force  vanishes  from  Eq.  (9)  if  Pen  =  Pex-l  As  in  the  previous  section,  we  simply 
drop  the  inertial  term  in  (9)  and  replace  Pen  by  niTe,  so  that  En  is  specified  by  the  equation 


-eE„  = 


B  a  ntTeii 
njaC  IBI 


-aB 

+  [t-^  +  VemeUe||. 


(10) 


The  determination  of  the  transverse  electric  field  Ex  involves  the  ion  dynamics  and,  in  the 
case  of  a  bounded  plasma  in  a  conducting  vessel,  also  couples  to  the  sheath  potentials. 
This  will  be  discussed  in  a  subsequent  publication. 

3.  Formal  Analysis  of  Mode  Structure  and  Stability 

A.  Linearized  Normal  Modes 

To  elucidate  the  way  in  which  the  electric  field  from  Eq.  (4)  or  Eq.  (10)  couples 
the  electron  and  ion  densities,  we  shall  examine  the  linear  normal  modes  supported  by  the 
system.  Assuming  an  equilibrium  with  uniform  density  no  and  temperature  To, 
linearizing  Eq.  (4),  assuming  normal  modes  of  the  form  and  neglecting  collisions, 

we  have 


-eE  =  ikTe 


+ 


ikTp 

no 


n: 


(11) 
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In  this  analysis,  we  use  linearized  cold  fluid  equations  for  the  ions,  so  that 


n;  = 


ikngeE 

co^mi 


(12) 


Combining  Eqs.  (11)  and  (12),  we  find 


eE  =  - 


ikT. 


i-kV/o" ' 


(13) 


To  complete  the  analysis,  we  use  the  linearized  Vlasov  equation  for  the  electrons. 


34  npeEaFp 

3t  9z  m.  dv 


(14) 


where  Fo(v)  is  the  normalized  equilibrium  electron  velocity  distribution  and  fe(z,v,t)  is 
now  the  first-order  perturbation.  Using  (13)  in  (14),  we  find 


T 

f.  =-— 


Fq  (v) 


m 


e  l-k%^/0)^  v-(o/k' 


(15) 


Te  can  be  calculated  by  using 


ne=/dvfe(v) 


(16) 


and 

neTo + noTe  =  J  dv  mv^  fe(v).  (17) 

Using  (16)  and  (17)  in  (15)  yields  a  dispersion  relation 


8 


(18) 


kV  f  dv(v^-Ve^)Fo  (v) 

J  v-ra/k 

We  note  first  that  if  Cs  ->0,  so  that  the  ions  are  immobile,  then 

(0  =  ±kVe  (^^) 

is  a  pair  of  exact  solutions  of  Eq.  (18).  Here,  Ve  =  /  m^  is  the  electron  thermal 

velocity.  This  mode  represents  the  unphysical  high  frequency  oscillations  that  keep  the 
electron  density  closely  coupled  to  the  ion  density.  However,  the  oscillation  frequency  o 
is  much  smaller  than  co^.  In  a  simulation  with  spatial  grid  scale  Az,  the  fiiequency  is 
limited  to  (B<  Ve/Az,  and  if  there  is  any  spatial  smoothing  the  highest-k  modes  are  strongly 
damped,  so  that  the  highest  meaningful  frequency  is  actually  much  less.  In  a  typical 
application  such  as  our  ECR  simulations,  Az  may  be  about  1  cm,  and  Ve/Az  of  the  order  of 
10*  sec'\  as  compared  to  (Oj,  of  order  lO"  sec'‘. 

We  next  examine  the  mode  structure  when  c*  *  0.  To  perform  the  integrals  simply 
in  closed  form,  we  assume  that  Fo(v)  is  a  Lorentzian  distribution, 


Fo(v)  =  7  2 1'  — 
Jtv  +v/ 


Then  Eq.  (18)  becomes 


(20) 


1- 


2 


dxx(l-x^) 
_^(l  +  x^)^(x-Q) 


(21) 
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where  x  s  v/ve  and  Q  =  (O  /kve.  The  contour  for  the  integral  in  Eq.  (21)  can  be  closed 
above,  and  according  to  causality  the  pole  at  x  =  £2  is  to  be  enclosed  in  the  contour,  even 
if  Im  Q  is  negative.  Using  the  method  of  residues,  Eq.  (21)  reduces  to 


i^J _ 2Q^-iQ(l-£2^) 

i  £2^  ~  (1  +  £2^)2 


If  we  multiply  out  the  denominators  in  (22),  we  obtain  a  sixth  order  poljmomial 
equation. 


£2^(1  +  £2^  )^  -  4£2^  +  2i£2^  (1  -  £2^ )  =  —  (1  +  £2^  )^ . 

ni: 


(23) 


We  know  that  two  of  the  roots  lie  at  £2  =  ±1  if  me/mi-»0.  Assuming  that  these  two  roots 
lie  close  to  ±1,  a  perturbative  solution  gives 


f 


CO  =  kv. 


±l  +  i 


V 


(24) 


These  modes  are  thus  seen  to  be  slightly  unstable,  by  order  me/mi.  This  is  not  significant; 
the  mode  growth  is  so  slow  that  it  is  obliterated  by  any  of  a  number  of  incoherence  effects 
(and  in  fact  the  entire  formalism  is  based  on  neglect  of  higher  order  in  me/mj.)  For 
example,  we  shall  see  in  the  next  subsection  that  even  the  slightest  degree  of  spatial 
smoothing  damps  the  mode. 

Next,  we  look  for  low  frequency  modes  with  l£2l«l.  Keeping  only  lowest  order 
in  £2  in  the  real  and  imaginary  terms  of  Eq.  (23),  we  find  one  pair  of  low-frequency  roots. 


(0  =  kc, 


( 

m. 

±l-iJ 

i 

mi 

) 

(25) 
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These  are  the  ion  sound  modes,  with  the  correct  dispersion  relation  in  the  limit  X.d  — >  0, 
and  with  the  correct  representation  of  electron  Landau  damping  for  the  Lorentzian 
distribution  (20).  We  show  in  the  Appendix  that  when  a  complete  Vlasov  treatment  is 
used,  both  electron  and  ion  Landau  damping  terms  appear  correctly,  and  ion  sound 
instability  is  also  correctly  represented  if  there  is  electron-ion  streaming. 

The  last  two  roots  of  the  sixth-order  equation  (23)  are  a  double  root  at  =  i. 
However,  this  is  a  spurious  root,  which  is  not  a  root  of  the  original  equation  (22).  Thus 
the  formalism  supports  only  two  pairs  of  modes.  One  pair  is  the  ion  sound  mode, 
correctly  represented.  The  second  pair  of  modes  are  the  (unphysical,  but  essentially 
stable)  high  frequency  modes  which  tightly  couple  n®  to  nj  and  thus  preserve 
quasineutrality. 

B.  Effect  of  Spatial  Smoothing 

hi  practice,  we  find  that  it  is  necessary  to  apply  some  spatial  smoothing  to  the 
electric  field,^’  to  overcome  the  fluctuations  introduced  by  particle  statistics  and  enhanced 
by  the  derivative  operation  in  Eq.  (4).  One  might  wonder  about  the  effect  of  smoothing 
on  the  stability  of  the  scheme.  Thus  we  reexamine  the  mode  structure  in  the  presence  of 
smoothing. 

The  smoothed  electric  field  E  may  be  represented  as 
E(z)  =  jr..dz'K(z-z')E(zO,  (26) 

where  E(z)  is  given  by 


^  8Te  To  ans 
-eE(z)  =  -Tr^+“ — - 


dz 


no  dz 


(27) 
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and  K(z-zO  is  some  symmetric  kernel  normalized  to  unity.  In  Fourier  representation,  the 
convolution  becomes  simply 


—  eE|j  =  — eEijK]j  =  ikKjjTg  +  ikKj^  nj . 


no 


(28) 


Note  that  the  Fourier  transform  Kk  is  always  real,  and  1  >  IKkl  for  all  k.  If  the  width  of  the 
smoothing  kernel  K(z)  is  narrow  compared  to  k'*,  then  Kk  -^1  and  the  effect  of  smoothing 
vanishes.  Using  (28)  in  place  of  (1 1),  we  can  retrace  the  derivation  of  the  dispersion 
relation  (22).  We  obtain  the  modified  form 


1-Kk 


1  _nir  -iQ(l-Q^) 

^  (1  +  a^)^ 


(29) 


Equation  (29)  can  be  solved  in  the  same  way  as  (22).  Again,  there  are  four  genuine  roots. 
The  two  high-frequency  modes  which  couple  the  electrons  to  the  ions  are  now 


0>k=kVe 


±  ^(2-Ky)Ky,  -  i(l  -  Kk)  +  order 


m: 


(30) 


These  modes  are  now  seen  to  be  damped  as  long  as  1-Kk  >  order  (me/mO,  so  the  electrons 
should  follow  the  ions  in  a  quiescent  fashion.  The  two  low-frequency  solutions  of  (29) 
are 


05k  =kc,Kk 


1/2 


±l-iKk 


(31) 


Thus,  as  might  be  expected,  the  ion  sound  waves  experience  a  reduction  in  frequency  and 
in  Landau  damping,  if  there  is  smoothing  on  a  scale  comparable  to  the  wavelength. 


12 


Obviously,  if  one  wishes  to  resolve  sound  waves  of  a  given  wavelength,  smoothing 
should  only  be  applied  on  smaller  spatial  scales. 

4.  Numerical  Implementation 

The  choice  of  time  steps  is  limited  by  a  number  of  considerations,  in  addition  to 
the  obvious  requirement  that  the  time  step  resolve  any  time  scale  of  interest,  such  as  a 
wave  period,  (i)  Accuracy  requires  that  during  a  single  time  step,  particles  not  traverse  a 
range  over  which  the  electric  field,  or  other  macroscopic  variables,  change  significantly, 
(ii)  If  collisions  are  an  important  aspect  of  the  problem,  and  Monte  Carlo  methods  are 
used  to  model  them,  the  time  steps  for  each  species  must  also  be  limited  to  a  fraction  of 
the  collision  time.  Both  of  these  conditions  typically  allow  the  ion  time  step  to  be  longer 
than  the  electron  time  step  by  a  factor  of  the  order  of  (mj/noe)*^.  (iii)In  addition,  we  have 
seen  that  Eq.  (4)  couples  the  electrons  to  the  ions  by  inducing  rapid  stable  electron 
oscillations  with  phase  velocity  equal  to  the  electron  thermal  velocity  Vc.  Since  these 
oscillations  can  be  excited  by  statistical  fluctuations  in  a  single  cell,  a  conservative 
procedure  to  avoid  numerical  difficulties  is  to  choose  the  electron  time  step  no  larger  than 
the  cell  size  divided  by  Ve,  and  to  recalculate  the  electric  field  acting  on  the  electrons  at 
each  electron  time  step,  hi  practice,  we  have  found  that  these  conditions  should  be 
satisfied  in  the  most  stressing  situations,  such  as  the  collisionless  simulations  presented 
later  in  this  paper.  Collisional  situations  generally  are  more  forgiving,  and  appear  often  to 
allow  longer  time  steps.  Thus,  to  summarize,  we  use  relatively  long  ion  time  steps, 
chosen  to  satisfy  conditions  (i)  and  (ii),  and  subdivide  these  time  steps,  typically  by  a 
factor  of  the  order  of  (mi/me)’^,  to  obtain  electron  time  steps  that  satisfy  all  three 
conditions.  These  conditions  permit  electron  time  steps  that  are  typically  three  orders  of 
magnitude  larger  than  the  time  steps  At  <  2/oOj,  that  are  needed  for  conventional  PIC  codes. 
The  ion  time  steps  are  even  larger,  and  in  addition  the  spatial  grid  scales  can  be  orders  of 
magnitude  larger  than  Xd- 

The  electric  field  is  calculated  as  a  grid  quantity  at  each  electron  time  step,  from 
Eq.  (4),  or  Eq.  (10)  if  the  plasma  is  magnetized.  This  electric  field  is  then  applied  directly 
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to  the  electrons  at  each  time  step.  To  push  the  ions  over  a  long  ion  time  step,  an  electric 
field  is  used  which  is  simply  the  average  of  the  electric  fields  at  each  of  the  electron  time 
steps  during  this  ion  time  step.  Thus  there  is  significant  temporal  smoothing  of  Eq.  (4)  or 
(10),  primarily  a  smoothing  of  the  fluctuations  in  T*,  which  helps  to  reduce  statistical 
fluctuations  in  the  ion  motion. 

The  electrostatic  potential  energy  of  the  plasma  is  jdV  (n|  -  n^)  ,  and  therefore 

is  zero  to  the  extent  that  exact  quasineutrality  is  maintained.  Thus,  in  the  absence  of 
inelastic  scattering,  the  total  kinetic  energy  should  be  conserved.  This  presumes,  of 
course,  that  the  electrons  and  ions  see  exactly  the  same  electric  field.  However,  energy 
conservation  is  preserved  even  if  the  electrons  are  pushed  subject  to  the  instantaneous 
electric  field,  and  the  ions  to  a  time-averaged  field,  as  long  as  the  time  averaging  is 
properly  centered  over  the  ion  time  step.  In  practice,  the  oscillation  of  the  electrons  about 
the  ion  density  profile  means  that  quasineutrality  is  not  exactly  satisfied  at  any  given  time 
step,  and  correlations  are  possible  which  lead  to  long-term  drifts  in  the  total  kinetic 
energy.  However,  these  are  very  slow  and  can  be  controlled  by  giving  proper  attention  to 
maintaining  a  sufficient  number  of  simulation  particles  and  adequate  spatial  smoothing. 
Some  examples  will  be  given. 

For  the  collisionless  examples  below,  we  have  used  a  large  number  of  particles 
per  grid  point,  ranging  from  an  average  of  500  particles  per  cell  for  free  expansion  to 
24(K)  per  cell  for  the  ion  acoustic  wave.  However,  in  our  2D  ECR  simulation  code,  where 
there  is  significant  electron  collisionality,  we  get  good  results  with  about  200  particles  per 
cell.  We  use  a  linear  laydown  for  both  the  electrons  and  the  ions.  The  ion  density  and  the 
pressure  are  smoothed  using  standard  filtering  techniques.*  The  algorithm  for  pushing 
the  particles  is  a  centered  difference  scheme.  The  electrons  are  subcycled  with  typically 
32  electron  time  steps  for  each  ion  time  step.  For  the  examples  shown  here,  the  system  is 
periodic,  but  in  our  ECR  code  the  same  method  is  used  for  a  bounded  system. 
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5.  An  Example:  Free  Expansion  of  a  Plasma 

To  illustrate  the  use  of  the  quasineutral  formulation,  we  consider  a  test  problem 
which  can  be  solved  exactly,  and  also  can  be  solved  analytically  within  our  formulation. 
Consider  the  free  expansion  of  a  plasma  consisting  of  hot  isothermal  electrons  and  cold 
ions,  beginning  with  a  Gaussian  spatial  profile  in  one  dimension: 


fi(z,v,0) 


8(v), 


(32a) 


fe(z,v,0) 


1  nie 


2^ 

meV 
2T()g  ^ 


(32b) 


Even  though  this  is  a  problem  that  is  easily  solved  analytically,  it  poses  a  stiff  challenge 
to  a  particle  simulation,  since  there  is  a  wide  range  in  plasma  density,  and  we  assume 
there  are  no  collisions.  (Collisions  make  it  much  easier  to  implement  this  type  of 
technique,  by  smoothing  out  statistical  fluctuations.) 


A.  Analytic  Treatment:  Exact  Quasineutrality 

An  exact  solution  to  the  Vlasov  equation,  with  initial  conditions  specified  by  Eqs. 
(32),  can  be  found  in  closed  form.  It  corresponds  to  self-similar  isothermal  expansion. 


fi(z,v,t) 


VjcLi(t) 


exp 


zL,(t)^ 
Li(t)  / 


(33a) 


fe(z.V,t) 


1  me 

ViLe(t)\27lTe(t) 


z"  m  f  zLe(t)^ 
Le^Ct)  2Te(t)l,''  Le(t)^ 


(33b) 
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The  expansion  is  driven  by  the  electron  pressure,  with  the  ions  dragged  along  by  the 
electrostatic  field.  Thus  a  complete  solution,  using  Poisson’s  equation,  would  show  that 
Le(t)  is  larger  than  Li(t)  by  a  small  amount  of  the  order  of  the  Debye  length.  But  to  the 
extent  that  quasineutrality  is  observed,  Le(t)  =  Li(t),  and  Li(t)  and  Te(t)  are  determined  by 
conservation  of  momentum  and  energy  as  the  solutions  to 


Li(t)  = 


2T,(t) 

miLi(t)  ’ 


(34a) 


Te(t)  =  T^(0)  - 1  miL;  2(t) .  (34b) 

Equations  (34)  can  be  solved  for  Li(t)  in  closed  form.  For  the  initial  conditions  Specified 
by  Eq.  (32),  the  solution  is 


Li(t)  =  VLo^+2c3oV. 


(34c) 


B.  Analytic  Treatment  Using  Our  Model 

Within  our  model,  with  E  given  by  Eq.  (4),  it  is  easy  to  show  that  the  electrons 
and  ions  each  expand  self-similarly,  as  in  Eqs.  (33),  but  with  Le(t)  not  exactiy  equal  to 
Li(t).  Li(t)  and  Tc  (t)  are  given  by  Eqs.  (34),  but  Le(t)  is  given  by 


Le=2Ve^ 


vLe 


h 

u 


(35) 


Rewriting  Eq.  (35)  in  terms  of  A  =  (Le-Li)/Li,  and  subtracting  (34a)  from  (35)  gives 


2VeV(2  +  A)A  ^  m,^ 
[  1  + A  mj/ 


(36) 
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The  first  term  of  Eq.  (36)  causes  A  to  oscillate  about  an  equilibrium  point.  The  second 
term  causes  a  very  small  offset  to  the  equilibrium  point.  (Interestingly,  this  offset  is 
negative,  so  that  in  this  formulation  the  ions  lead  the  electrons  by  a  displacement  of  order 
nie/mi.)  Thus  it  is  reasonable  to  assume  IAI«  1  and  simplify  (36)  to 


A 


^2Ve(t)'>" 

L|(t) 


A. 


(37) 


Equation  (37)  indicates  that  A  oscillates  at  the  rapid  frequency  (Oo=2ve/Li(t),  and  since  coo 
is  fast  compared  to  the  time  scale  for  ion  motion,  these  can  be  considered  to  be  simple 
harmonic  oscillations  at  a  slowly  varying  frequency.  These  oscillations  are  not  physical; 
they  ate  the  mechanism  for  coupling  the  electrons  to  the  ions  within  our  quasineutral 
formulation  (4).  However,  they  are  stable  oscillations  which  normally  maintain  a  very 
small  amplitude  (in  fact,  comparable  to  the  statistical  fluctuations  that  would  otherwise  be 
present  due  to  the  finite  number  of  particles),  and  thus  are  of  no  real  significance.  As 
mentioned  earlier,  a  conservative  condition  for  numerical  stability  is  that  the  time  step  be 
less  than  coo'S  and  in  fact,  since  statistical  fluctuation  can  occur  on  length  scales  down  to 
a  single  cell  size  Az,  less  than  Az/Vc. 

C.  Numerical  Simulation 

Figures  1-3  show  the  results  of  a  numerical  simulation  of  the  free  expansion 
problem,  for  a  hydrogen  plasma  (mi/ntic  =  1836).  However,  it  was  not  possible  in  the 
simulation  to  allow  the  plasma  to  expand  into  a  true  vacuum,  since  the  electric  field  from 
Eq.  (4)  depends  on  the  reciprocal  of  the  plasma  density.  (The  noisy  electric  field  in 
regions  with  very  low  density  eventually  dominates  the  problem.)  In  order  to  control  the 
noise,  we  added  a  low  density  floor  (5%  of  the  peak  density)  to  the  Gaussian  density 
profile  as  can  be  seen  in  Fig.  la.  We  use  200,000  macroparticle  electrons  and  an  equal 
number  of  ions,  in  a  ID  system  400  cells  long..  The  cell  size  is  Az  =  1  cm.  The  system  is 
initialized  in  accordance  with  Eq.  (32),  with  Lo  =  14  cm  and  Te  =  0.33  eV.  The  time  step 
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is  At  =  25  ns.  Figure  1  shows  the  electron  density  profile  (solid  curve)  and  the  ion  density 
profile  (dashed  curve)  at  t  =  0  and  t  =  20  ps.  The  deviations  from  quasi-neutrality  are 
visible  only  at  the  peak  density,  and  at  the  boundary  with  the  low  density  background 
plasma.  We  note  that  the  expansion  is  very  nearly  self-similar,  as  predicted.  (Eventually, 
small  deviations  from  self-similarity,  due  to  the  presence  of  the  low-density  background, 
become  visible.)  Figure  2  shows  Li(t)  from  the  simulation.  The  x’s  show  the  analytic 
solution  (34c)  for  Li(t).  We  see  that  the  expansion  is  smooth,  and  the  analytic  solution  is 
well  verified.  Figure  3  shows  plots  of  the  total  electron  kinetic  energy  We(t),  the  total  ion 
kinetic  energy  Wi(t),  and  the  total  kinetic  energy  W(t).  We  note  that  the  effect  of  the 
expansion  is  to  convert  the  electron  kinetic  energy  (which  is  nearly  all  thermal)  to  ion 
streaming  energy.  Overall  kinetic  energy  is  conserved  to  within  4%.  To  the  extent  that 
quasineutrality  is  maintained,  there  is  essentially  no  potential  energy,  since  the  electron 
potential  energy  is  always  exactly  the  negative  of  the  ion  potential  energy. 


6.  Second  Example:  Ion  Sound 

Ion  sound  with  wavelength  much  greater  than  Xd  is  a  simple  example  of  a  plasma 
mode  that  occurs  on  the  ion  time  scale  and  maintains  quasineutrality.  It  is  therefore  not 
an  easy  mode  to  simulate  with  a  standard  PIC  code  using  Poisson’s  equation  and  a 
realistic  value  of  nie/mi.  In  Sec.  3  above,  we  used  a  simplified  model  (Vlasov  electrons 
with  Lorentzian  distribution,  cold  fluid  ions)  to  show  that  the  ion  sound  mode  is 
contained  within  our  quasineutral  model,  and  in  the  Appendix  we  provide  a  full  Vlasov 
analysis  that  shows  that  the  model  gives  exactly  the  correct  dispersion  relation,  including 
electron  and  ion  Landau  damping  terms.  Here  we  show  the  results  of  simulations  of  ion 
sound  within  the  nearly-linear  and  nonlinear  regimes,  with  the  mass  ratio  mi/me  =  1836. 
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A.  Standlng'Wave  Ion  Sound  in  tiie  Linear  Be^me 

We  initiate  an  i<Mi  sound  standing  wave  by  loading  initial  particle  densities 

nesnisi)o[l-f  asin(2icz/L)], 

with  the  fluid  velocities  Ue  and  U|  initially  zero  for  both  species.  The  simulation  is  done 
with  periodic  boundary  conditions  in  a  ID  ^tem  of  length  2L= 25  cm,  with  cell  length 
Az  ss  0.125  cm;  thus,  there  are  two  wavelengths  witlun  the  box,  and  100  cells  per 
wavelength.  The  wave  amplitude  is  taken  to  be  a  *  0.05,  and  2400  particles  of  each 
species  are  used,  per  cell.  A  very  large  number  of  particles  is  needed  for  this  simulation 
to  resolve  tiie  very  small  wave  amplitude.  The  electron  teiiq>erature  is  set  to  133  eV,  so 
that  Ve=4.8xl0^  cto/scc,  Cg=  l.lx  10^  cm/scc,  and  the  wave  period  should  be  Uct=  113 
jis.  The  ions  are  initially  cold.  The  ion  time  step  is  15.6  ns,  and  the  electrons  are 
subcycled  4  times  per  ion  time  step. 

The  temporal  evolution  of  the  fundamental  (wavelength  L)  Fourier  mode  of  rie  and 
m  is  shown  in  Fig.  4.  The  densities  are  taken  from  the  Fourier  transforms  of  the 
instantaneous  values  of  ne(z)  and  ni(z)  at  the  ion  time  step.  We  see  that  the  wave  period  is 
11.1  |is,  in  excellent  agreement  with  the  theoretical  value.  We  note  that  at  any  given  time, 
quasineutrality  (n*  =  nO  holds  to  better  than  5%.  Plots  of  the  spatial  dependence  of  the 
wave  (not  shown)  indicate  that  the  two  wavelengths  contained  within  the  box  are  exactly 
the  same. 

Equation  (A9)  predicts  Landau  damping  at  the  rate  9%  per  wave  period. 
However,  in  a  one-dimensional  ion  sound  wave,  even  at  wave  anqrlitude  o(  =  0.05, 
trapping  of  the  resonant  electrons  interferes  with  Landau  damping  at  a  very  early  stage. 
The  trapping  fiequency  is  *=  V^sTm^,  and  according  to  the  linearized  version  of 
Eq.  (4),  E  =  kofTe .  Thus 

©X  *=  kCj^a—  =  9.6  kCj ,  OT 
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and  Landau  damping  should  be  over  in  a  fraction  of  a  wave  period.  The  slight  damping 
seen  in  Fig.  4  thus  would  seem  to  be  the  initial  Landau  damping,  followed  by  a  slower 
damping  due  to  mode  coupling  to  higher  harmonics.  The  second  harmonic  (wavelength 
L/2)  grows  to  a  level  of  10%  of  the  fundamental  after  two  periods,  as  shown  in  Fig.  4. 
This  is  in  reasonable  agreement  with  simple  two-mode  coupling  theory,  which  predicts 
15%  at  this  time,  and  consequent  damping  of  the  fundamental  by  about  1.5%.  From  a 
computational  point  of  view,  it  is  worth  noting  that  the  rate  of  mode  coupling  is  found  to 
be  quite  sensitive  to  the  grid  spacing.  The  simulation  of  Fig.  4,  with  100  grid  points  per 
wavelength,  shows  qualitative  agreement  with  theory,  but  a  calculation  with  25  grid 
points  per  wavelength  shows  substantially  faster  second  harmonic  growth  and  decay  of 
the  fundamental. 

B.  Nonlinear  li-aveling  Wave 

The  nonlinear  evolution  of  a  standing  wave  is  complex,  but  for  traveling  sound 

waves  an  analytic  solution  is  possible  via  the  method  of  characteristics.  Since  Ad  is  set 

to  zero  within  our  model,  the  ion  sound  waves  are  nondispersive,  and  therefore  the  mode 
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coupling  theory  is  essentially  the  same  as  that  of  ordinary  sound  waves  in  a  neutral  gas, 
with  Te  playing  the  role  of  the  gas  temperature.  (However,  ion  sound  waves  are 
isothermal,^^  whereas  sound  waves  in  a  molecular  gas  are  adiabatic.  Therefore  the 
adiabatic  constant  y  must  be  set  to  unity  in  the  ion  sound  theory.)  The  velocity  is  found  to 
be  the  solution  of  the  implicit  equation 

u(z,  t)  =  u(z  -  u(z,  t)) .  (40) 

It  is  well  known  that  Eq.  (40)  leads  to  steepening  of  the  density  and  velocity  profiles,  and 
ultimately  to  wave  breaking  when  8u/9z  becomes  infinite.  For  an  initial  sinusoidal 
profile 
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u(z,0) 


=  asm 


(41) 


this  occurs  when 


L 

t-tbreak-2jjojc^- 


(42) 


In  Fig.  5  we  show  the  results  of  a  simulation  similar  to  those  of  the  previous  subsection, 
except  that  the  wave  amplitude  is  larger  at  a  =  0.20  and  a  traveling  wave  is  initiated  by 
setting  the  initial  ion  fluid  velocity  as  in  Eq.  (41).  We  note  the  steepening  of  the  wave  up 
to  the  point  of  breaking  at  about  t  =  9  psec,  in  good  agreement  with  the  theoretical 
prediction  t  =  9.04  psec.  As  the  wave  nears  the  breaking  point,  non-physical  structure 
such  as  the  ripple  and  sharp  peak  in  Fig.  5b  begin  to  appear.  These  features  appear  to  be 
associated  with  finite  spatial  resolution,  and  limit  the  accuracy  of  determination  of  the 
breaking  point  by  perhaps  1  psec.  After  wave  breaking,  the  quasineutral  theory  is  no 
longer  correct,  since  A-d  length  scales  are  relevant  to  the  subsequent  evolution  of  the  ion 
acoustic  shocks. 


7.  Conclusions 

We  have  presented  a  method  for  doing  plasma  simulations  with  particle  electrons 
and  particle  ions,  in  the  quasineutral  limit.  The  method  permits  (indeed,  it  requires)  the 
use  of  grid  spacing  long  compared  to  A.d  and  time  stepping  long  compared  to  We 
have  demonstrated  analytically  that  the  method  is  stable  (at  least  on  a  continuous  time 
basis),  and  that  it  gives  the  correct  dispersion  relation  for  ion  sound,  including  Landau 
resonance  terms.  We  have  also  demonstrated  the  use  of  the  method  to  simulate  free 
expansion  of  a  plasma,  and  linear  and  nonlinear  ion  sound. 
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We  believe  that  the  technique  presented  here  should  be  useful  for  a  variety  of 
problems  involving  quasineutral  plasmas.  We  are  using  this  method  in  a  2D  simulation 
code  for  magnetized  plasmas,  which  includes  in  addition  cross-field  transport,  sheaths  at 
insulating  or  conducting  walls,  ECR  heating,  and  collisions  and  chemistry.  The 
techniques  used  to  model  these  other  aspects  of  the  physics  will  be  discussed  in 
subsequent  publications.  The  code  runs  for  very  long  times  (hundreds  of  (isec),  with  time 
steps  typically  on  the  order  of  10  nsec,  and  with  running  time  typically  on  the  order  of  a 
few  hours  on  a  IBM  RS6000  workstation. 


22 


Appendix:  Vlasov  Theory  of  Ion  Sound 


In  this  section  we  show  analytically  that  our  quasineutral  formulation,  using  the 
approximate  form  (4)  for  E,  gives  exactly  the  correct  dispersion  relation  for  linearized  ion 
sound  waves,  including  even  the  Landau  damping  terms. 

Review;  Vlasov-Poisson  Theory  of  Ion  Sound 

To  clarify  the  issues  involved,  we  begin  with  a  brief  review  of  the  standard 
derivation  of  the  ion  sound  dispersion  relation  from  the  linearized  ID  Vlasov  equation, 

=  0,  (Al) 

0t  dz  m„  3v 

with  the  upper  sign  for  electrons  and  the  lower  sign  for  ions,  and  with  the  electric  field  E 
determined  by  Poisson’s  equation.  Using  a  linearized  normal  mode  representation  for  the 
perturbed  part  of  the  distribution  function  fa,  we  find 


fak(v)  =  ± 


inneE  dF^  1 
m„  dv  (D-kv  ’ 


(A2) 


and  the  species  densities  no*  are 


dvFoa  (v) 
m-kv 


Inserting  (A3)  into  Poisson’s  equation. 


nc*  =  ± 


inpeE 

m„ 


ikEk  =  47c(nek-nik), 


(A3) 


(A4) 


we  obtain  the  dispersion  relation. 
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0  =  1- 


(A5) 


Fo/(v) 

v-co/k 


Foi'(v) 

V  — ©/  k 


where  causality  indicates  that  the  contours  in  Eqs.  (A3)  and  (A5)  go  below  the  pole  at  v  = 
co/k.  The  usual  dispersion  relation  for  ion  sound  can  be  obtained  by  making  several 
approximations  in  Eq.  (A5).  First  assume  that  the  imaginary  part  of  o/k  is  small.  Let  u 
(which  may  be  zero)  be  the  relative  electron-ion  drift,  and  assume  that  Foi(v)  is  a 
synunetric  function  of  v,  but  Foe  is  a  symmetric  function  of  w  =  v  -  u.  Assume  that  Te » 
T|,  u  «  Ve,  and  Vi «  Ci)/k«  Vc.  Then  the  first  integral  in  Eq.  (A5)  can  be  treated  by  first 
extracting  the  contribution  from  the  pole  at  v  =  oo/k,  and  then  expanding  the  remaining 
non-singular  integral  under  the  assumption  that  code  «  Iwl  and  u  «  Iwl  for  nearly  all 
electrons: 


The  second  integral  in  Eq.  (5)  can  be  treated  by  first  extracting  the  contribution  from  the 
pole,  and  then  expanding  the  remaining  integral  under  the  assumption  that  code »  Ivl  for 
nearly  all  ions: 


Keeping  only  the  leading  term  in  the  expansion  (A6)  gives  the  dispersion  relation. 
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0  =  1- 


k^ca 


p' 


CO 


■fdw^Foi 


+  (0 


(A8) 

In  the  case  of  Maxwellian  velocity  distributions,  this  reduces  to  the  familiar  form 


kCj 


Mt 


kCc 


2  d  +  k^Xo") 


2x3/2 


T  ' 

—  Fn, 


L^e 


ifle 


(A9) 


Quasineutral  Derivation 

In  the  quasineutral  context,  the  electric  field  is  determined  by  the  linearized  form 
ofEq.(4), 


On  Un 


(AlO) 


with 


Pek  =  J  dvfek(v)m^v^  =  inocE 


dvFo,  (V) 

(D-kv 


(All) 


Using  Eqs.  (A3)  and  (All)  in  (AlO),  we  obtain  the  dispersion  relation 


f  dv(v^  -  Ve^) 

IF 

/ 

Oe  2 

dvFoi' 

J 

v-(o/k 

—  Cg 

v-co/k 

(A12) 


Equation  (A12)  looks  different  from  Eq.  (A8).  Nevertheless,  using  the  expansions  (A6) 
and  (A7)  but  keeping  all  terms  shown  in  (A6)  gives 
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0  =  - 


CO„i 

'"pi 


k^co 


p* 


0 


Jdw’Fo,-^ 


(D 


pc 


+  (0 


P> 


(A13) 


Equation  (A13)  is  identical  to  Eq.  (A8),  except  for  the  absence  of  the  first  term  on  the 
RHS  of  (A8),  This  term  is  smaller  than  the  other  terms  in  (A8)  by  order  and  thus 

disappears  in  the  quasineutral  limit  kXo  0. 
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width  (cm) 


time  (|X  sec) 

Fig.  2  —  Simulation  result  for  expansion  of  the  ion  characteristic  with  Li(t). 
The  x’s  are  the  analytic  result  (34c). 
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Fig.  3  —  Time  dependence  of  the  electron  kinetic  energy  W,(t)  (dot-dashed  curve),  ion  kinetic  energy 
Wi(t)  (dashed  curve),  and  total  kinetic  energy  W(t)  (solid  curve)  during  free  expansion. 
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Fig.  5  —  Density  profiles  showing  steepening  of  the  traveling  ion  sound  wave: 
(a)  t  =  0,  (b)  t  =  5.75  /^sec,  (c)  t  =  9.0  jitsec. 
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Fig.  5  (Continued)  —  Density  profiles  showing  steepening  of  the  traveling  ion  sound  wave: 
(a)  t  =  0,  (b)  t  =  5.75  jiisec,  (c)  t  =  9.0  /isec. 
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Fig.  5  (Continued)  —  Density  profiles  showing  steepening  of  the  traveling  ion  sound  wave 
(a)  t  =  0,  (b)  t  =  5.75  fisec,  (c)  t  =  9.0  fxsec. 
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